home *** CD-ROM | disk | FTP | other *** search
/ Power Programmierung / Power-Programmierung (Tewi)(1994).iso / magazine / drdobbs / 1989 / 09 / mclaughl.lst < prev    next >
File List  |  1989-07-27  |  13KB  |  496 lines
  1. _SIMULATED ANNEALING_
  2. by Michael P. McLaughlin
  3.  
  4.  
  5. [LISTING ONE]
  6.  
  7. /* ANNEAL.C -- Author: Dr. Michael P. McLaughlin
  8. */
  9. #include <stdio.h>
  10. #include <math.h>
  11.  
  12. #define ABS(x) ((x<0)?(-(x)):(x))
  13. #define BOOLEAN int
  14. #define TRUE 1
  15. #define FALSE 0
  16. #define FORWARD 1
  17. #define BACK 0
  18.  
  19. struct cardtype {
  20.    char card[3];
  21.    int code;
  22.    } cards[25];
  23. float ratios[10][2];
  24. int tableau[25],best_tableau[25];
  25. float temperature,best_temperature,t_low,t_min,t_ratio;
  26. long seed,score,best_score,tot_successes,tot_failures,scor,exg_cutoff,
  27.      report_time,report_interval,modulus = 2147483399,step = 1;
  28. int next;
  29. BOOLEAN equilibrium,frozen = FALSE,quench=FALSE,final=FALSE;
  30.  
  31. /* variables needed to assess equilibrium */
  32. float totscore,totscore2,avg_score,sigma,half_sigma,score_limit;
  33. long bscore,wscore,max_change,nsuccesses,nfailures,scale,
  34.      incount,outcount;
  35. int incount_limit,outcount_limit,succ_min,first_limit,ultimate_limit;
  36.  
  37. /* poker hands in tableau */
  38. int hand[12][5] = {0,1,2,3,4,      /* rows */
  39.                    5,6,7,8,9,
  40.                    10,11,12,13,14,
  41.                    15,16,17,18,19,
  42.                    20,21,22,23,24,
  43.                    0,5,10,15,20,   /* columns */
  44.                    1,6,11,16,21,
  45.                    2,7,12,17,22,
  46.                    3,8,13,18,23,
  47.                    4,9,14,19,24,
  48.                    0,6,12,18,24,   /* diagonals */
  49.                    4,8,12,16,20};
  50.                 /* 0,4,12,20,24     = corners */
  51.  
  52. long rand1()
  53. /* Get uniform 31-bit integer -- Ref. CACM, June 1988, pg. 742 */
  54. {
  55.    register long k;
  56.    k = seed/52774;
  57.    seed = 40692*(seed-k*52774)-k*3791;
  58.    if (seed < 0)
  59.       seed += modulus;
  60.    return(seed);
  61. }
  62. double uni(z)
  63. int z;
  64. {
  65.    return((double) z*rand1()/modulus);
  66. }
  67. void initialize()
  68. /* Set up entire simulated annealing run. */
  69. {
  70.    char label[20];
  71.    double exgc;
  72.    register i;
  73.    FILE *fp;
  74.    fp = fopen("a.dat","r");
  75.    fscanf(fp,"%f %s\n",&temperature,label);
  76.    fscanf(fp,"%f %s\n",&t_low,label);
  77.    fscanf(fp,"%f %s\n",&t_min,label);
  78.    fscanf(fp,"%f %s\n",&avg_score,label);
  79.    fscanf(fp,"%ld %s\n",&scale,label);
  80.    fscanf(fp,"%f %s\n",&t_ratio,label);
  81.    fscanf(fp,"%f %s\n",&sigma,label);
  82.    fscanf(fp,"%lf %s\n",&exgc,label);
  83.    fscanf(fp,"%d %s\n",&succ_min,label);
  84.    fscanf(fp,"%d %s\n",&incount_limit,label);
  85.    fscanf(fp,"%d %s\n",&outcount_limit,label);
  86.    fscanf(fp,"%d %s\n",&first_limit,label);
  87.    fscanf(fp,"%d %s\n",&ultimate_limit,label);
  88.    fscanf(fp,"%ld %s\n",&report_interval,label);
  89.    fscanf(fp,"%ld %s\n",&seed,label);
  90.    half_sigma = 0.5*sigma; report_time = report_interval;
  91.    score_limit = 20; exg_cutoff = modulus*exgc;
  92.    for (i=0;i<10;i++)
  93.       fscanf(fp,"%f %f\n",&ratios[i][0],&ratios[i][1]);
  94.    for (i=0;i<25;i++) {
  95.       fscanf(fp,"%d %s\n",&cards[i].code,cards[i].card);
  96.       tableau[i] = cards[i].code;
  97.    }
  98.    fclose(fp);
  99.    printf("                   TOTAL     TOTAL  CURRENT AVERAGE\n");
  100.    printf("STEP TEMPERATURE SUCCESSES FAILURES  SCORE   SCORE ");
  101.    printf(" SIGMA\n\n");
  102. }
  103. void init()
  104. /* set up for new temperature */
  105. {
  106.    nsuccesses = 0; nfailures = 0; equilibrium = FALSE; incount = 0;
  107.    outcount = 0; bscore = 0; wscore = 100000; max_change = 0;
  108.    totscore = 0; totscore2 = 0;
  109.    if (temperature < t_min)
  110.       frozen = TRUE;
  111. }
  112. void get_new_temperature()
  113. {
  114.    if (temperature < ratios[next][0]) {
  115.       t_ratio = ratios[next][1];
  116.       next++;
  117.    }
  118.    temperature = t_ratio*temperature;
  119. }
  120. void check_equilibrium()
  121. /* Determine whether equilibrium has been reached. */
  122. {
  123.    if ((nsuccesses+nfailures) > ultimate_limit)
  124.       equilibrium = TRUE;
  125.    else if (nsuccesses >= succ_min) {
  126.       if (incount > incount_limit)
  127.          equilibrium = TRUE;
  128.       else {
  129.          if (outcount > outcount_limit) {
  130.             if (nsuccesses > first_limit)
  131.                equilibrium = TRUE;
  132.             else {
  133.                incount = 0;
  134.                outcount = 0;
  135.             }
  136.          }
  137.       }
  138.    }
  139. }
  140. void update()
  141. /* Compute statistics, etc. */
  142. {
  143.    float ascore,s;
  144.    if (nsuccesses > 0) {
  145.       ascore = totscore/nsuccesses;
  146.       avg_score = ascore+scale;
  147.       s = totscore2/nsuccesses-ascore*ascore;
  148.       if (s > 0.0) {
  149.          sigma = sqrt(s); half_sigma = 0.5*sigma;
  150.          score_limit = avg_score-half_sigma;
  151.       }
  152.    }
  153.    if ((temperature < t_low)&&((bscore-wscore) == max_change))
  154.       frozen = TRUE;
  155. }
  156. void selectwr(array,mode,nchoices)
  157. /* Select from array without replacement. */
  158. int *array,mode,nchoices;
  159. {
  160.    int i,temp,size,index;
  161.    size = mode==1 ? 12 : 25;
  162.    for (i=0;i<nchoices;i++) {
  163.       index = uni(size);
  164.       temp = array[--size];
  165.       array[size] = array[index];
  166.       array[index] = temp;
  167.    }
  168. }
  169. void reconfigure(direction)
  170. /* If direction is FORWARD, make a new move; if it is BACK, restore
  171.    the tableau to its previous configuration. */
  172. int direction;
  173. {
  174.    static long partition[6] = {0,0,715827799,1296031795,1766307855,
  175.                                2147483400};
  176.    static int hindex[12] = {0,1,2,3,4,5,6,7,8,9,10,11};
  177.    static int cindex[25] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,
  178.                             14,15,16,17,18,19,20,21,22,23,24};
  179.    static int overlap[66] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,0,5,10,15,
  180.                              20,1,6,11,16,21,-1,2,7,12,17,22,-1,-1,3,
  181.                              8,13,18,23,-1,-1,-1,4,9,14,19,24,-1,-1,-1,
  182.                              -1,0,6,12,18,24,0,6,12,18,24,4,8,12,16,20,
  183.                              20,16,12,8,4,12};
  184.    static int mode,nchoices,common,last;   /* remember changes */
  185.    int temp,hi,lo,i,j,k;
  186.    if (direction == FORWARD) {
  187.       if (rand1() < exg_cutoff) {   /* giant step */
  188.          mode = 1; nchoices = 2;
  189.          selectwr(hindex,mode,nchoices);
  190.          hi = hindex[11]; lo = hindex[10];   /* order hand indices */
  191.          if (hi < lo) {
  192.             hi = lo;
  193.             lo = hindex[11];
  194.          }
  195.          common = overlap[hi*(hi-1)/2+lo];   /* triangular matrix */
  196.       }
  197.       else {   /* baby step */
  198.          mode = 0; nchoices = 2;
  199.          rand1();   /* How many cards to rotate? */
  200.          while (seed > partition[nchoices])
  201.             nchoices++;
  202.          selectwr(cindex,mode,nchoices);
  203.          temp = tableau[cindex[24]];   /* rotate forward */
  204.          for (i=1,j=24;i<nchoices;i++,j--)
  205.             tableau[cindex[j]] = tableau[cindex[j-1]];
  206.          tableau[cindex[j]] = temp;
  207.          last = j;
  208.       }
  209.    }
  210.    if (mode == 1) {   /* swapping hands is symmetrical */
  211.       register first,second,c;
  212.       first = hindex[11]; second = hindex[10];
  213.       k = (common<0) ? 5 : 4;
  214.       for (c=0,i=0,j=0;c<k;c++,i++,j++) {
  215.          if (hand[first][i] == common)
  216.             i++;
  217.          if (hand[second][j] == common)
  218.             j++;
  219.          temp = tableau[hand[first][i]];
  220.          tableau[hand[first][i]] = tableau[hand[second][j]];
  221.          tableau[hand[second][j]] = temp;
  222.       }
  223.    }
  224.    else if (direction == BACK) {   /* rotation is not */
  225.       temp = tableau[cindex[last]];
  226.       for (i=1,j=last;i<nchoices;i++,j++)
  227.          tableau[cindex[j]] = tableau[cindex[j+1]];
  228.       tableau[cindex[24]] = temp;
  229.    }
  230. }
  231. long get_score(crds)
  232. /* Return score of one hand. */
  233. int *crds;
  234. {
  235.    static long scores[18] = {0,1,1,20,1,9,20,1760,1,9,9,
  236.                              293,20,293,1760,108,215,27456};
  237.    int flush,i,index = 0;
  238.    flush = crds[0]&crds[1]&crds[2]&crds[3]&crds[4]&0xff00;
  239.    for (i=0;i<5;i++)
  240.       crds[i] = crds[i]&0xff;   /* ignore suits */
  241.    {
  242.       register b,k,t;   /* sort cards from high to low */
  243.       for (b=4;b>0;b--)
  244.          for (k=0;k<b;k++)
  245.             if (crds[k] < crds[k+1]) {
  246.                t = crds[k];
  247.                crds[k] = crds[k+1];
  248.                crds[k+1] = t;
  249.             }
  250.    }
  251.    if (!flush) {   /* multiple cards ? */
  252.       if (crds[0] == crds[1]) index += 8;
  253.       if (crds[1] == crds[2]) index += 4;
  254.       if (crds[2] == crds[3]) index += 2;
  255.       if (crds[3] == crds[4]) index += 1;
  256.    }
  257.    if (!index) {   /* straight ? */
  258.       if (((crds[0]-crds[4]) == 4)||((crds[0] == 14)&&(crds[1] == 5)))
  259.          index = 15;
  260.       if (flush)
  261.          index = (index == 15) ? 17 : 16;
  262.    }
  263.    return(scores[index]);
  264. }
  265. long evaluate1()
  266. /* Return total score of tableau (no corners). */
  267. {
  268.    int temp[5],h,c;
  269.    long scr = 0;
  270.    for (h=0;h<12;h++) {   /* h<13 for corners */
  271.       for (c=0;c<5;c++)
  272.          temp[c] = tableau[hand[h][c]];
  273.       scr += get_score(temp);
  274.    }
  275.    return(scr);
  276. }
  277. void decide()
  278. /* Either accept new configuration or restore old configuration */
  279. {
  280.    float pdt1,pdt2,sscor;
  281.    long s;
  282.    BOOLEAN acceptable = FALSE;
  283.    scor = evaluate1();
  284.    if (scor >= score) {
  285.       if (scor > best_score) {
  286.          register i;
  287.          best_score = scor;
  288.          best_temperature = temperature;
  289.          for (i=0;i<25;i++)
  290.             best_tableau[i] = tableau[i];
  291.       }
  292.       acceptable = TRUE;
  293.    }
  294.    else if (temperature > 0.0) {   /* uphill movement? */
  295.       if (uni(1) < exp((scor-score)/temperature))   /* Equation 2 */
  296.          acceptable = TRUE;
  297.    }
  298.    if (acceptable) {   /* statistics, etc. */
  299.       s = ABS(score-scor);
  300.       if (s > max_change)
  301.          max_change = s;
  302.       score = scor;
  303.       sscor = scor-scale;   /* to aid precision of sigma */
  304.       totscore += sscor;
  305.       totscore2 += sscor*sscor;
  306.       nsuccesses++;
  307.       if (ABS(avg_score-score) < half_sigma)
  308.          incount++;
  309.       else
  310.          outcount++;
  311.       if (score > bscore)   /* maximization */
  312.          bscore = score;
  313.       else if (score < wscore)
  314.          wscore = score;
  315.    }
  316.    else {   /* unacceptable */
  317.       reconfigure(BACK);
  318.       nfailures++;
  319.    }
  320. }
  321. void report()
  322. {
  323.    tot_successes += nsuccesses; tot_failures += nfailures;
  324.    printf("%3ld %10.1f %9ld %9ld %7ld %7.0f  %5.1f\n",step,temperature,
  325.           tot_successes,tot_failures,score,avg_score,sigma);
  326.    report_time += report_interval;
  327.    if (frozen) {
  328.       int temp[25],k,kk;
  329.       register i,j;
  330.       BOOLEAN ok;
  331.       if (final)
  332.          printf("\nFINAL -- ");
  333.       else if (quench)
  334.          printf("\nQUENCH -- ");
  335.       else
  336.          printf("\nFROZEN -- ");
  337.       printf("BEST SCORE = %5ld  BEST TEMPERATURE = %5.1f\n\n",
  338.               best_score,best_temperature);
  339.       for (i=0;i<25;i++) {
  340.          k = best_tableau[i];
  341.          j = 0;
  342.          ok = FALSE;
  343.          while (!ok) {
  344.             if (cards[j].code == k) {
  345.                temp[i] = j;
  346.                ok = TRUE;
  347.             }
  348.             else j++;
  349.          }
  350.       }
  351.       for (i=0;i<25;i+=5) {
  352.          for (j=i;j<(i+5);j++) {
  353.             k = 4-strlen(cards[temp[j]].card);
  354.             for (kk=0;kk<k;kk++) printf(" ");
  355.             printf("%s",cards[temp[j]].card);
  356.          }
  357.          printf("\n");
  358.       }
  359.       printf("\n");
  360.    }
  361. }
  362. void try_new_temperature()
  363. {
  364.    init();
  365.    while (!equilibrium) {
  366.       reconfigure(FORWARD);
  367.       decide();
  368.       check_equilibrium();
  369.    }
  370.    update();
  371.    if (!quench) {   /* ratchet */
  372.       while ((float) score < score_limit) {   /* maximization */
  373.          reconfigure(FORWARD);
  374.          decide();
  375.       }
  376.    }
  377. }
  378. main()
  379. {
  380.    initialize();
  381.    while (!frozen) {
  382.       try_new_temperature();
  383.       update();
  384.       if ((step == report_time)||(frozen))
  385.          report();
  386.       step++;
  387.       if (temperature > 0.0)
  388.          get_new_temperature();
  389.    }
  390. /* Quench frozen configuration. */
  391.    temperature = 0;
  392.    quench = TRUE;
  393.    try_new_temperature();
  394.    update();
  395.    report();
  396.    step++;
  397. /* Quench best configuration (rarely useful). */
  398.    {
  399.       register i;
  400.       final = TRUE;
  401.       for (i=0;i<25;i++)
  402.          tableau[i] = best_tableau[i];
  403.    }
  404.    try_new_temperature();
  405.    update();
  406.    report();
  407.    printf("SEED = %ld\n",seed);   /* seed of next run */
  408. }
  409.  
  410.  
  411.  
  412. [LISTING TWO]
  413.  
  414.      if ((nsuccesses+nfailures) > ULTIMATE_LIMIT)
  415.        equilibrium = true;
  416.      else if (nsuccesses >= SUCC_MIN) {
  417.        if (incount > INCOUNT_LIMIT)
  418.         equilibrium = true;
  419.        else {
  420.         if (outcount > OUTCOUNT_LIMIT) {
  421.          if (nsuccesses > FIRST_LIMIT)
  422.           equilibrium = true;
  423.          else {
  424.           incount = 0;
  425.           outcount = 0;
  426.          }
  427.         }
  428.        }
  429.      }
  430.  
  431.  
  432. [LISTING THREE]
  433.  
  434. 2150 temperature
  435. 1.5 t_low
  436. 0.1 t_min
  437. 40 avg_score
  438. 4400 scale
  439. 0.7 t_ratio
  440. 150 sigma
  441. 0.2 exg_cutoff
  442. 200 succ_min
  443. 250 incount_limit
  444. 400 outcount_limit
  445. 2700 first_limit
  446. 10000 ultimate_limit
  447. 1 report_interval
  448. 1234567890 seed
  449. 360 0.8
  450. 215 0.7
  451. 85 0.8
  452. 60 0.9
  453. 30 0.95
  454. 15 0.9
  455. 7 0.8
  456. 3 0.7
  457. 0 0.7
  458. 0 0.7
  459. 1032 8H
  460. 2061 KS
  461. 270 AC
  462. 526 AD
  463. 522 10D
  464. 1029 5H
  465. 265 9C
  466. 1035 JH
  467. 1037 KH
  468. 525 KD
  469. 515 3D
  470. 263 7C
  471. 2058 10S
  472. 2062 AS
  473. 518 6D
  474. 1036 QH
  475. 267 JC
  476. 259 3C
  477. 2053 5S
  478. 269 KC
  479. 520 8D
  480. 1028 4H
  481. 1038 AH
  482. 519 7D
  483. 1026 2H
  484.  
  485.  
  486.  
  487. [EXAMPLE 1] 
  488.  
  489. Equation 1:  N /su/hi//N /su/lo/ = exp((E lo-E hi)/kT)     
  490.  
  491. Equation 2:  Prob(accept worse score) = exp((S /su/lo/ -S /su/hi/)/T)   
  492.  
  493.  
  494. Example 1: Equations that describe the Boltzman distribution   
  495.  
  496.